C $TITLE: 'FBNGF'
C $NOFLOATCALLS
      SUBROUTINE FBNGF (NEQ,NEQ2,IRESRV,IB11,IC11,ID11,IX11,IW)
C     FBNGF SETS THE BLOCKING PARAMETERS FOR THE B, C, AND D ARRAYS
C     FOR OUT-OF-CORE STORAGE.
      INTEGER*4 NEQ,NEQ2,IB11,IC11,ID11,IX11,NBLN,NDLN,NBCD,IR,IRESX
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
C**
C     D      WRITE(*,*) '  FBNGF: START'
C**
      IRESX=IRESRV-IMAT
      NBLN=NEQ*NEQ2
      NDLN=NEQ2*NEQ2
      NBCD=2*NBLN+NDLN
      IF (NBCD.GT.IRESX) GO TO 1
      ICASX=1
      IB11=IMAT+1
      GO TO 2
1      CONTINUE
C**
C     D      WRITE(*,*) '  FBNGF: OPEN 11,12,13,14,15'
C**
      OPEN (11,FORM='UNFORMATTED')
      OPEN (12,FORM='UNFORMATTED')
      OPEN (13,FORM='UNFORMATTED')
      OPEN (14,FORM='UNFORMATTED')
      OPEN (15,FORM='UNFORMATTED')
      IF (ICASE.LT.3) GO TO 3
      IF (NBCD.GT.IRESRV.OR.NBLN.GT.IRESX) GO TO 3
      ICASX=2
      IB11=1
2     NBBX=1
      NPBX=NEQ
      NLBX=NEQ
      NBBL=1
      NPBL=NEQ2
      NLBL=NEQ2
      GO TO 5
3     IR=IRESRV
      IF (ICASE.LT.3) IR=IRESX
      ICASX=3
C**
C     D      WRITE(*,*) '  FBNGF: OPEN 16'
C**
      OPEN (16,FORM='UNFORMATTED')
      IF (NDLN.GT.IR) ICASX=4
      NBCD=2*NEQ+NEQ2
      NPBL=IR/NBCD
      NLBL=IR/(2*NEQ2)
      IF (NLBL.LT.NPBL) NPBL=NLBL
      IF (ICASE.LT.3) GO TO 4
      NLBL=IRESX/NEQ
      IF (NLBL.LT.NPBL) NPBL=NLBL
4     IF (NPBL.LT.1) GO TO 6
      NBBL=(NEQ2-1)/NPBL
      NLBL=NEQ2-NBBL*NPBL
      NBBL=NBBL+1
      NBLN=NEQ*NPBL
      IR=IR-NBLN
      NPBX=IR/NEQ2
      IF (NPBX.GT.NEQ) NPBX=NEQ
      NBBX=(NEQ-1)/NPBX
      NLBX=NEQ-NBBX*NPBX
      NBBX=NBBX+1
      IB11=1
      IF (ICASE.LT.3) IB11=IMAT+1
5     IC11=IB11+NBLN
      ID11=IC11+NBLN
      IX11=IMAT+1
      WRITE(IW,11)  NEQ2
      IF (ICASX.EQ.1) RETURN
      WRITE(IW,8)  ICASX
      WRITE(IW,9)  NBBX,NPBX,NLBX
      WRITE(IW,10)  NBBL,NPBL,NLBL
C**
C     D      WRITE(*,*) '  FBNGF: RETURN AT END'
C**
      RETURN
6     WRITE(IW,7)  IRESRV,IMAT,NEQ,NEQ2
      WRITE(*,7)  IRESRV,IMAT,NEQ,NEQ2
      STOP
C
7     FORMAT(' ERROR - INSUFFICIENT STORAGE FOR INTERACTION MATRICIES',
     1'  IRESRV,IMAT,NEQ,NEQ2 =',4I5)
8     FORMAT(' FILE STORAGE FOR NEW MATRIX SECTIONS -  ICASX =',I2)
9     FORMAT(' B FILLED BY ROWS-',15X,'NO. BLOCKS =',I3,3X,
     1'ROWS PER BLOCK =',I3,3X,'ROWS IN LAST BLOCK =',I3)
10    FORMAT(' B BY COLUMNS, C AND D BY ROWS -',2X,'NO. BLOCKS =',
     1 I3,4X,'R/C PER BLOCK =',I3,4X,'R/C IN LAST BLOCK =',I3)
11    FORMAT (//,' N.G.F. - NUMBER OF NEW UNKNOWNS IS',I4)
      END
C
C
C
      SUBROUTINE FBLOCK (NROW,NCOL,IMAX,IRNGF,IPSYM,IW)
C     FBLOCK SETS PARAMETERS FOR OUT-OF-CORE SOLUTION FOR THE PRIMARY
C     MATRIX (A)
      INTEGER*4 IMX1,NROW,NCOL,IPSYM,KK,I,J,K,KA,NOP
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      REAL*8 ARG,PHAZ
      COMPLEX*16 SSX,DETER
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,
     1 IMAT,ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON /SMAT/ SSX(16,16)
C     D      WRITE(*,*) '  FBLOCK: START'
      IMX1=IMAX-IRNGF
      IF (NROW*NCOL.GT.IMX1) GO TO 2
      NBLOKS=1
      NPBLK=NROW
      NLAST=NROW
      IMAT=NROW*NCOL
      IF (NROW.NE.NCOL) GO TO 1
      ICASE=1
C     D      WRITE(*,*) '  FBLOCK: RETURN BEFORE 1'
      RETURN
1     ICASE=2
      GO TO 5
2      CONTINUE
C**
C     D      WRITE(*,*) '  FBLOCK: OPEN 11,12,13,14'
C**
      OPEN (11,FORM='UNFORMATTED')
      OPEN (12,FORM='UNFORMATTED')
      OPEN (13,FORM='UNFORMATTED')
      OPEN (14,FORM='UNFORMATTED')
      IF (NROW.NE.NCOL) GO TO 3
      ICASE=3
      NPBLK=IMAX/(2*NCOL)
      NPSYM=IMX1/NCOL
      IF (NPSYM.LT.NPBLK) NPBLK=NPSYM
      IF (NPBLK.LT.1) GO TO 12
      NBLOKS=(NROW-1)/NPBLK
      NLAST=NROW-NBLOKS*NPBLK
      NBLOKS=NBLOKS+1
      NBLSYM=NBLOKS
      NPSYM=NPBLK
      NLSYM=NLAST
      IMAT=NPBLK*NCOL
      WRITE(IW,14)  NBLOKS,NPBLK,NLAST
      GO TO 11
3     NPBLK=IMAX/NCOL
      IF (NPBLK.LT.1) GO TO 12
      IF (NPBLK.GT.NROW) NPBLK=NROW
      NBLOKS=(NROW-1)/NPBLK
      NLAST=NROW-NBLOKS*NPBLK
      NBLOKS=NBLOKS+1
      WRITE(IW,14)  NBLOKS,NPBLK,NLAST
      IF (NROW*NROW.GT.IMX1) GO TO 4
      ICASE=4
      NBLSYM=1
      NPSYM=NROW
      NLSYM=NROW
      IMAT=NROW*NROW
      WRITE(IW,15)
      GO TO 5
4     ICASE=5
      NPSYM=IMAX/(2*NROW)
      NBLSYM=IMX1/NROW
      IF (NBLSYM.LT.NPSYM) NPSYM=NBLSYM
      IF (NPSYM.LT.1) GO TO 12
      NBLSYM=(NROW-1)/NPSYM
      NLSYM=NROW-NBLSYM*NPSYM
      NBLSYM=NBLSYM+1
      WRITE(IW,16)  NBLSYM,NPSYM,NLSYM
      IMAT=NPSYM*NROW
5      CONTINUE
      NOP=NCOL/NROW
      IF (NOP*NROW.NE.NCOL) GO TO 13
      IF (IPSYM.GT.0) GO TO 7
C
C     SET UP SSX MATRIX FOR ROTATIONAL SYMMETRY.
C
      PHAZ=6.2831853072D0/NOP
      DO 6 I=2,NOP
      DO 6 J=I,NOP
      ARG=PHAZ*(I-1)*(J-1)
      SSX(I,J)=DCMPLX(DCOS(ARG),DSIN(ARG))
      SSX(J,I)=SSX(I,J)
6      CONTINUE
      GO TO 11
C
C     SET UP SSX MATRIX FOR PLANE SYMMETRY
C
7     KK=1
      SSX(1,1)=DCMPLX(1.,0.)
      IF ((NOP.EQ.2).OR.(NOP.EQ.4).OR.(NOP.EQ.8)) GO TO 8
      STOP
8     KA=NOP/2
      IF (NOP.EQ.8) KA=3
      DO 10 K=1,KA
      DO 9 I=1,KK
      DO 9 J=1,KK
      DETER=SSX(I,J)
      SSX(I,J+KK)=DETER
      SSX(I+KK,J+KK)=-DETER
9     SSX(I+KK,J)=DETER
10    KK=KK*2
11      CONTINUE
C**
C     D      WRITE(*,*) '  FBLOCK: RETURN AFTER 11'
C**
      RETURN
12    WRITE(IW,17)  NROW,NCOL
      WRITE(*,17)  NROW,NCOL
      STOP
13    WRITE(IW,18)  NROW,NCOL
      WRITE(*,18)  NROW,NCOL
      STOP
C
14    FORMAT (//' MATRIX FILE STORAGE-NO. BLOCKS=',I5,
     1' COLUMNS PER BLOCK=',I5,' COLUMNS IN LAST BLOCK=',I5)
15    FORMAT(' SUBMATRICIES FIT IN CORE')
16    FORMAT(' SUBMATRIX PARTITIONING-NO. BLOCKS=',I5,
     1' COLUMNS PER BLOCK=',I5,' COLUMNS IN LAST BLOCK=',I5)
17    FORMAT(' ERROR - INSUFFICIENT STORAGE FOR MATRIX',2I5)
18    FORMAT(' SYMMETRY ERROR - NROW,NCOL=',2I5)
      END
C
C
C
      SUBROUTINE LOAD(ZARRAY,ZLR,ZLI,ZLC,SI,BI,LD,ITAG,LDTYP,
     1 LDTAG,LDTAGF,LDTAGT,IW)
C
C     LOAD CALCULATES THE IMPEDANCE OF SPECIFIED SEGMENTS FOR VARIOUS
C     TYPES OF LOADING
C
C***
C***  FIX TO ADD PERMEABILITY TO WIRES - 12 DEC 88 - R W ADLER
C***  ADD VALUE OF REL. PERM. AFTER CONDUCTIVITY ON LD CARD
C***  2 LINES ADDED - 1 LINE CHANGED - ALSO FIXES TO FUNCTION ZINT
C***
      CHARACTER ADUM*20
      REAL*8 TPCJX,ZREAL,ZIMAG
      COMPLEX*16 TPCJ
      COMPLEX*16 ZARRAY,ZT,ZINT
      INTEGER*4 ITAG,N1,N2,N,NP,M1,M2,M,MP,IPSYM
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      COMMON /ZLOAD/ NLOAD,NLODF
      DIMENSION LDTYP(30),LDTAG(30),LDTAGF(30),LDTAGT(30),ZLR(30),
     1 ZLI(30),ZLC(30),TPCJX(2),ZARRAY(LD),ITAG(LD),SI(LD),BI(LD)
      EQUIVALENCE (TPCJ,TPCJX)
      DATA TPCJX/0.,1.883698955D+9/
C**
C     D      WRITE(*,*) '  LOAD: START, NLOAD=',NLOAD,' NLODF=',NLODF
      FZERO=0.
C**
C
C     WRITE(IW,HEADING)
C
      WRITE(IW,25)
C
C     INITIALIZE D ARRAY, USED FOR TEMPORARY STORAGE OF LOADING
C     INFORMATION.
C
      DO 1 I=N2,N
      ZARRAY(I)=DCMPLX(0.,0.)
1      CONTINUE
      IWARN=0
C
C     CYCLE OVER LOADING CARDS
C
      ISTEP=0
 2    ISTEP=ISTEP+1
      IF (ISTEP.LE.NLOAD) GO TO 5
      IF (IWARN.EQ.1) WRITE(IW,26)
      IF (N1+2*M1.GT.0) GO TO 4
      NOP=N/NP
      IF (NOP.EQ.1) GO TO 4
      DO 3 I=1,NP
      ZT=ZARRAY(I)
      L1=I
      DO 3 L2=2,NOP
      L1=L1+NP
      ZARRAY(L1)=ZT
3      CONTINUE
 4    CONTINUE
C**
C     D      WRITE(*,*) '  LOAD: RETURN AFTER 4'
C**
      RETURN
 5    IF (LDTYP(ISTEP).LE.5) GO TO 6
      WRITE(IW,27)  LDTYP(ISTEP)
      STOP
 6    LDTAGS=LDTAG(ISTEP)
      JUMP=LDTYP(ISTEP)+1
      ICHK=0
C
C     SEARCH SEGMENTS FOR PROPER ITAGS
C
      L1=N2
      L2=N
      IF (LDTAGS.NE.0) GO TO 7
      IF (LDTAGF(ISTEP).EQ.0.AND.LDTAGT(ISTEP).EQ.0) GO TO 7
      L1=LDTAGF(ISTEP)
      L2=LDTAGT(ISTEP)
      IF (L1.GT.N1) GO TO 7
      WRITE(IW,29)
      STOP
7      CONTINUE
      DO 17 I=L1,L2
      IF (LDTAGS.EQ.0) GO TO 8
      IF (LDTAGS.NE.ITAG(I)) GO TO 17
      IF (LDTAGF(ISTEP).EQ.0) GO TO 8
      ICHK=ICHK+1
      IF (ICHK.GE.LDTAGF(ISTEP).AND.ICHK.LE.LDTAGT(ISTEP)) GO TO 9
      GO TO 17
 8    ICHK=1
C
C     CALCULATION OF LAMDA*IMPED. PER UNIT LENGTH, JUMP TO APPROPRIATE
C     SECTION FOR LOADING TYPE
C
9      CONTINUE
C**
C     D      WRITE(*,*) '  LOAD: AFTER  9, JUMP=',JUMP
C**
      GO TO (10,11,12,13,14,15), JUMP
10      CONTINUE
      ZT=ZLR(ISTEP)/SI(I)+TPCJ*ZLI(ISTEP)/(SI(I)*WLAM)
      IF(ABS(ZLC(ISTEP)).GT.1.E-20) ZT=ZT+WLAM/(TPCJ*SI(I)*ZLC(ISTEP))
      GO TO 16
 11   ZT=TPCJ*SI(I)*ZLC(ISTEP)/WLAM
      IF(ABS(ZLI(ISTEP)).GT.1.E-20) ZT=ZT+SI(I)*WLAM/(TPCJ*ZLI(ISTEP))
      IF (ABS(ZLR(ISTEP)).GT.1.E-20) ZT=ZT+SI(I)/ZLR(ISTEP)
      ZT=1./ZT
      GO TO 16
 12   ZT=ZLR(ISTEP)*WLAM+TPCJ*ZLI(ISTEP)
      IF (ABS(ZLC(ISTEP)).GT.1.E-20) ZT=ZT+1./(TPCJ*SI(I)*SI(I)*
     1 ZLC(ISTEP))
      GO TO 16
13      ZT=TPCJ*SI(I)*SI(I)*ZLC(ISTEP)
      IF (ABS(ZLI(ISTEP)).GT.1.E-20) ZT=ZT+1./(TPCJ*ZLI(ISTEP))
      IF (ABS(ZLR(ISTEP)).GT.1.E-20) ZT=ZT+1./(ZLR(ISTEP)*WLAM)
      ZT=1./ZT
      GO TO 16
14      ZT=CMPLX(ZLR(ISTEP),ZLI(ISTEP))/SI(I)
      GO TO 16
C***
C***  FIX FOR PERM. OF WIRES
C***
C15     CONTINUE
C15      RMU = ZLI(STEP)
15      RMU = ZLI(ISTEP)
        IF(RMU.EQ.0)RMU = 1.
C       ZT=ZINT(ZLR(ISTEP)*WLAM,BI(I))
        ZT=ZINT(ZLR(ISTEP)*WLAM,RMU,BI(I))
C***
C***  END FIX
C***
16      CONTINUE
      ZREAL=ABS(DREAL(ZARRAY(I)))
      ZIMAG=ABS(DIMAG(ZARRAY(I)))
      IF ((ZREAL+ZIMAG).GT.1.D-20) IWARN=1
      ZARRAY(I)=ZARRAY(I)+ZT
17      CONTINUE
      IF (ICHK.NE.0) GO TO 18
      WRITE(IW,28)  LDTAGS
      STOP
C
C     PRINTING THE SEGMENT LOADING DATA, JUMP TO PROPER PRINT
C
18      GOTO(19,20,21,22,23,24), JUMP
19      CONTINUE
      ADUM=' SERIES             '
      IDUM=2
      CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GO TO 2
20      CONTINUE
      ADUM='PARALLEL            '
      IDUM=2
      CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GOTO 2
21      CONTINUE
      ADUM='SERIES (PER METER)  '
      IDUM=5
      CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GOTO 2
22      CONTINUE
      ADUM='PARALLEL (PER METER)'
      IDUM=5
      CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GOTO 2
23      CONTINUE
      ADUM='FIXED IMPEDANCE     '
      IDUM=4
      CALL PRNT(FZERO,FZERO,FZERO,ZLR(ISTEP),ZLI(ISTEP),FZERO,
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GOTO 2
24      CONTINUE
      ADUM='  WIRE              '
      IDUM=2
      CALL PRNT(FZERO,FZERO,FZERO,FZERO,FZERO,ZLR(ISTEP),
     1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
      GOTO 2
C
25      FORMAT(//,7X,'LOCATION',10X,'RESISTANCE',3X,'INDUCTANCE',2X,
     1 'CAPACITANCE',7X,'IMPEDANCE (OHMS)',5X,'CONDUCTIVITY',4X,
     2 'TYPE',/,4X,'ITAG',' FROM THRU',10X,'OHMS',8X,'HENRYS',7X,
     3 'FARADS',8X,'REAL',6X,'IMAGINARY',4X,'MHOS/METER')
26      FORMAT(/,10X,'NOTE, SOME OF THE ABOVE SEGMENTS HAVE BEEN LOADED',
     1' TWICE - IMPEDANCES ADDED')
27      FORMAT(/,10X,'IMPROPER LOAD TYPE CHOOSEN, REQUESTED TYPE IS ',
     1 I3)
28      FORMAT(/,10X,'LOADING DATA CARD ERROR, NO SEGMENT HAS AN ITAG ',
     1 I5)
29      FORMAT(' ERROR - LOADING MAY NOT BE ADDED TO SEGMENTS IN N.G.F.',
     1 ' SECTION')
      END
C
C
C
      SUBROUTINE PRNT(FL1,FL2,FL3,FL4,FL5,FL6,IN1,IN2,IN3,ICHAR,
     1 IW,IA)
C
C     PRNT SETS UP THE PRINT FORMATS FOR IMPEDANCE LOADING
C
      CHARACTER ADUM*20,HALL*4,IFORM*6,IVAR*6,FMT*78,IA*4
      INTEGER*4 INTT
      DIMENSION IVAR(13),IA(5),IFORM(8),IN(3),INTT(3),FL(6),FLT(6)
      EQUIVALENCE (FMT,IVAR(1))
      DATA IFORM/'(/3X,','I5,','5X,','A5,','E13.4,','13X,','3X,',
     1 '5A4)'/
C
C     NUMBER OF CHARACTERS PER REAL*4 VARIABLE IS 4
C
      DATA HALL/' ALL'/
C**
C     D      WRITE(*,*) '   PRNT: START'
C**
      IN(1)=IN1
      IN(2)=IN2
      IN(3)=IN3
      FL(1)=FL1
      FL(2)=FL2
      FL(3)=FL3
      FL(4)=FL4
      FL(5)=FL5
      FL(6)=FL6
C
C     INTEGER FORMAT
C
      NINT=0
      IVAR(1)=IFORM(1)
      K=1
      I1=1
      IF (.NOT.(IN1.EQ.0.AND.IN2.EQ.0.AND.IN3.EQ.0)) GO TO 1
      INTT(1)=HALL
      NINT=1
      I1=2
      K=K+1
      IVAR(K)=IFORM(4)
1     DO 3 I=I1,3
      K=K+1
      IF (IN(I).EQ.0) GO TO 2
      NINT=NINT+1
      INTT(NINT)=IN(I)
      IVAR(K)=IFORM(2)
      GO TO 3
2     IVAR(K)=IFORM(3)
3     CONTINUE
      K=K+1
      IVAR(K)=IFORM(7)
C
C     FLOATING POINT FORMAT
C
      NFLT=0
      DO 5 I=1,6
      K=K+1
      IF (ABS(FL(I)).LT.1.E-20) GO TO 4
      NFLT=NFLT+1
      FLT(NFLT)=FL(I)
      IVAR(K)=IFORM(5)
      GO TO 5
4     IVAR(K)=IFORM(6)
5     CONTINUE
      K=K+1
      IVAR(K)=IFORM(7)
      K=K+1
      IVAR(K)=IFORM(8)
      WRITE(IW,FMT) (INTT(I),I=1,NINT),(FLT(J),J=1,NFLT),(IA(L),L=1,
     1 ICHAR)
C**
C     D      WRITE(*,*) '   PRNT: RETURN  IW=',IW
C**
      RETURN
      END
C
C
C
      SUBROUTINE GFOUT(CM,ZARRAY,X,Y,Z,SI,BI,ALP,BET,SALP,
     1 ICON1,ICON2,ITAG,IP,IW,IGFL,LD,LD2,IRESRV)
C
C     WRITE N.G.F. FILE TO 'TAPE.[IGFL]'
C
CLARGE: CM
      COMPLEX CM
      COMPLEX*16 ZARRAY,SSX
      COMPLEX*16 ZRATI,ZRATI2,T1,FRATI
C**
      REAL*4 DXA,DYA,XSA,YSA
      COMPLEX*8 AR1,AR2,AR3,EPSCF
C**
      INTEGER*4 ICON1,ICON2,ITAG,N1,N2,N,NP,M1,M2,M,MP,IPSYM,
     1 IOUT,IDM1,NEQ,NOP,NPEQ
      INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
      COMMON /GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
     1 IFAR,IPERF,T1,T2
      COMMON /GGRID/AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),
     1 DYA(3),XSA(3),YSA(3),NXA(3),NYA(3)
      COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
     1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
      COMMON /SMAT/ SSX(16,16)
      COMMON /ZLOAD/ NLOAD,NLODF
      COMMON /SAVE/ KCOM,COM(19,5),EPSR,SIG,SCRWLT,SCRWRT,FMHZ
      DIMENSION CM(IRESRV),ZARRAY(LD)
      DIMENSION X(LD),Y(LD),Z(LD),SI(LD),BI(LD),ALP(LD),BET(LD),
     1 SALP(LD),ICON1(LD),ICON2(LD),ITAG(LD),IP(LD2)
C**      DATA IGFL/20/
C**
C $NODEBUG
C**
C     D      WRITE(*,*) '   GFOUT: START'
C**
      NEQ=N+2*M
C**
$DEBUG
C**
      NPEQ=NP+2*MP
      NOP=NEQ/NPEQ
      IDM1=1
50      CONTINUE
C**      OPEN(IGFL,FILE='TAPE.20',FORM='UNFORMATTED',STATUS='UNKNOWN',
C**  1 ERR=100)
100      CONTINUE
      WRITE(*,'(A,I2,A)') ' OPEN UNIT ',IGFL,' FOR N.G.F. OUTPUT FILE'
200      CONTINUE
C**
      WRITE (IGFL) N,NP,M,MP,WLAM,FMHZ,IPSYM,KSYMP,IPERF,NRADL,EPSR,
     1SIG,SCRWLT,SCRWRT,NLOAD,KCOM
      IF (N.EQ.0) GO TO 1
      WRITE (IGFL) (X(I),I=1,N),(Y(I),I=1,N),(Z(I),I=1,N)
      WRITE (IGFL) (SI(I),I=1,N),(BI(I),I=1,N),(ALP(I),I=1,N)
      WRITE (IGFL) (BET(I),I=1,N),(SALP(I),I=1,N)
      WRITE (IGFL) (ICON1(I),I=1,N),(ICON2(I),I=1,N)
      WRITE (IGFL) (ITAG(I),I=1,N)
      IF (NLOAD.GT.0) WRITE (IGFL) (ZARRAY(I),I=1,N)
1     IF (M.EQ.0) GO TO 2
      J=LD-M+1
      WRITE (IGFL) (X(I),I=J,LD),(Y(I),I=J,LD),(Z(I),I=J,LD)
      WRITE (IGFL) (SI(I),I=J,LD),(BI(I),I=J,LD),(ALP(I),I=J,LD)
      WRITE (IGFL) (BET(I),I=J,LD),(SALP(I),I=J,LD)
      WRITE (IGFL) (ICON1(I),I=J,LD),(ICON2(I),I=J,LD)
      WRITE (IGFL) (ITAG(I),I=J,LD)
2     WRITE (IGFL) ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT
      IF (IPERF.EQ.2) WRITE (IGFL) AR1,AR2,AR3,EPSCF,DXA,DYA,XSA,YSA,NXA
     1,NYA
      IF (NOP.GT.1) WRITE (IGFL) ((SSX(I,J),I=1,NOP),J=1,NOP)
      WRITE (IGFL) (IP(I),I=1,NEQ),COM
      IF (ICASE.GT.2) GO TO 3
      IOUT=NEQ*NPEQ
      WRITE (IGFL) (CM(I),I=1,IOUT)
      GO TO 12
3     IF (ICASE.NE.4) GO TO 5
      REWIND 13
      I=NPEQ*NPEQ
      DO 4 K=1,NOP
      READ (13) (CM(J),J=1,I)
4     WRITE (IGFL) (CM(J),J=1,I)
      REWIND 13
      GO TO 12
5     REWIND 13
      REWIND 14
      IF (ICASE.EQ.5) GO TO 8
      IOUT=NPBLK*NEQ*2
      DO 6 I=1,NBLOKS
      CALL BLCKIN (CM,IDM1,IOUT,1,201,13)
6     CALL BLCKOT (CM,IDM1,IOUT,1,202,IGFL)
      DO 7 I=1,NBLOKS
      CALL BLCKIN (CM,IDM1,IOUT,1,203,14)
7     CALL BLCKOT (CM,IDM1,IOUT,1,204,IGFL)
      GO TO 12
8     IOUT=NPSYM*NPEQ*2
      DO 11 IOP=1,NOP
      DO 9 I=1,NBLSYM
      CALL BLCKIN (CM,IDM1,IOUT,1,205,13)
9     CALL BLCKOT (CM,IDM1,IOUT,1,206,IGFL)
      DO 10 I=1,NBLSYM
      CALL BLCKIN (CM,IDM1,IOUT,1,207,14)
10    CALL BLCKOT (CM,IDM1,IOUT,1,208,IGFL)
11    CONTINUE
      REWIND 13
      REWIND 14
12    REWIND IGFL
      WRITE(IW,13)  IGFL,IMAT
      WRITE(* ,13)  IGFL,IMAT
C**
      CLOSE(IGFL)
C**
C     D      WRITE(*,*) '   GFOUT: CLOSE N.G.F. OUTPUT FILE AND RETURN'
C**
      RETURN
C
13    FORMAT(///,' **** NUMERICAL GREEN''S FUNCTION FILE ON TAPE.',I2,
     1' ****',/,8X,'MATRIX STORAGE -',I7,' COMPLEX NUMBERS',///)
      END
C
C
C
      SUBROUTINE COUPLE (IW,CUR,WLAM,LD,LD3,ITAG)
C
C     COUPLE COMPUTES THE MAXIMUM COUPLING BETWEEN PAIRS OF SEGMENTS.
C
      INTEGER*4 ITAG
      REAL*8 DB10,GMAX,DBC,C
CLARGE: CUR
      COMPLEX CUR
      COMPLEX*16 Y11A,Y12A,Y11,Y12,Y22,YL,YIN,ZL,ZIN,RHO
      COMPLEX*16 VQD,VSANT,VQDS
      COMMON /YPARM/ NCOUP,ICOUP,NCTAG(5),NCSEG(5),Y11A(5),Y12A(20)
      COMMON /VSORC/ VQD(30),VSANT(30),VQDS(30),IVQD(30),ISANT(30),IQDS(
     130),NVQD,NSANT,NQDS
      DIMENSION CUR(LD3),ITAG(LD)
C**
C     D      WRITE(*,*) '  COUPLE: START'
C**
      IF (NSANT.NE.1.OR.NVQD.NE.0) RETURN
      J=ISEGNO(NCTAG(ICOUP+1),NCSEG(ICOUP+1),LD,ITAG)
      IF (J.NE.ISANT(1)) RETURN
      ICOUP=ICOUP+1
      ZIN=VSANT(1)
      Y11A(ICOUP)=CUR(J)*WLAM/ZIN
      L1=(ICOUP-1)*(NCOUP-1)
      DO 1 I=1,NCOUP
      IF (I.EQ.ICOUP) GO TO 1
      K=ISEGNO(NCTAG(I),NCSEG(I),LD,ITAG)
      L1=L1+1
      Y12A(L1)=CUR(K)*WLAM/ZIN
1     CONTINUE
      IF (ICOUP.LT.NCOUP) RETURN
      WRITE(IW,6)
      NPM1=NCOUP-1
      DO 5 I=1,NPM1
      ITT1=NCTAG(I)
      ITS1=NCSEG(I)
      ISG1=ISEGNO(ITT1,ITS1,LD,ITAG)
      L1=I+1
      DO 5 J=L1,NCOUP
      ITT2=NCTAG(J)
      ITS2=NCSEG(J)
      ISG2=ISEGNO(ITT2,ITS2,LD,ITAG)
      J1=J+(I-1)*NPM1-1
      J2=I+(J-1)*NPM1
      Y11=Y11A(I)
      Y22=Y11A(J)
      Y12=.5*(Y12A(J1)+Y12A(J2))
      YIN=Y12*Y12
C      DBC=CABS(YIN)
      DBC=ZABS(YIN)
      C=DBC/(2.*DREAL(Y11)*DREAL(Y22)-DREAL(YIN))
      IF (C.LT.0..OR.C.GT.1.) GO TO 4
      IF (C.LT..01) GO TO 2
      GMAX=(1.-DSQRT(1.-C*C))/C
      GO TO 3
2     GMAX=.5*(C+.25*C*C*C)
3     RHO=GMAX*DCONJG(YIN)/DBC
      YL=((1.-RHO)/(1.+RHO)+1.)*DREAL(Y22)-Y22
      ZL=1./YL
      YIN=Y11-YIN/(Y22+YL)
      ZIN=1./YIN
      DBC=DB10(GMAX)
      WRITE(IW,7)  ITT1,ITS1,ISG1,ITT2,ITS2,ISG2,DBC,ZL,ZIN
      GO TO 5
4     WRITE(IW,8)  ITT1,ITS1,ISG1,ITT2,ITS2,ISG2,C
5     CONTINUE
C**
C     D      WRITE(*,*) '  COUPLE: RETURN AT END'
C**
      RETURN
C
6     FORMAT (///,36X,'- - - ISOLATION DATA - - -',//,6X,
     1'- - COUPLING BETWEEN - -',8X,'MAXIMUM',15X,'- - - FOR MAXIMUM',
     2' COUPLING - - -',/,12X,'SEG.',14X,'SEG.',3X,'COUPLING',4X,
     3'LOAD IMPEDANCE (2ND SEG.)',7X,'INPUT IMPEDANCE',/,2X,'TAG/SEG.',
     43X,'NO.',4X,'TAG/SEG.',3X,'NO.',6X,'(DB)',8X,'REAL',9X,'IMAG.',
     5 9X,'REAL',9X,'IMAG.')
7     FORMAT (2(1X,I4,1X,I4,1X,I5,2X),F9.3,2X,1P,2(2X,E12.5,1X,E12.5))
8     FORMAT (2(1X,I4,1X,I4,1X,I5,2X),'**ERROR** COUPLING IS NOT ',
     1'BETWEEN 0 AND 1. (=',1P,E12.5,1H))
      END
C
C
C
C***
C***  FIX TO ADD PERMEABILITY TO WIRES - 12 DEC 88 - R W ADLER
C***  3 LINES CHANGED
C***
C     FUNCTION ZINT(SIGL,ROLAM)
      FUNCTION ZINT(SIGL,RMU,ROLAM)
C
C     ZINT COMPUTES THE INTERNAL IMPEDANCE OF A CIRCULAR WIRE
C
C
      REAL*8 PI,POT,TP,TPCMW,CNX,BER,BEI
      COMPLEX*16 CN,ZINT
      COMPLEX FJ
      COMPLEX*16 CC1,CC2,CC3,CC4,CC5,CC6,CC7,CC8,CC9,CC10,CC11,CC12,
     1 CC13,CC14,TH,PH,F,G,BR1,BR2
      DIMENSION FJX(2), CNX(2), CCN(28)
      EQUIVALENCE (FJ,FJX),(CN,CNX)
      EQUIVALENCE (CC1,CCN(1)),(CC2,CCN(3)),(CC3,C
     1CN(5)),(CC4,CCN(7)),(CC5,CCN(9)),(CC6,CCN(11)),(CC7,CCN(13)),
     2(CC8,CCN(15)),(CC9,CCN(17)),(CC10,CCN(19)),(CC11,CCN(21)),(CC1
     32,CCN(23)),(CC13,CCN(25)),(CC14,CCN(27))
      DATA PI,POT,TP,TPCMU/3.1415926D0,1.5707963D0,6.2831853D0,
     1 2.368705D+3/
      DATA CMOTP/60.00/,FJX/0.,1./,CNX/.70710678D0,.70710678D0/
      DATA CCN/6.E-7,1.9E-6,-3.4E-6,5.1E-6,-2.52E-5,0.,-9.06E-5,-9.01E-5
     1,0.,-9.765E-4,.0110486,-.0110485,0.,-.3926991,1.6E-6,-3.2E-6,1.17E
     2-5,-2.4E-6,3.46E-5,3.38E-5,5.E-7,2.452E-4,-1.3813E-3,1.3811E-3,-6.
     325001E-2,-1.E-7,.7071068,.7071068/
      TH(D)=(((((CC1*D+CC2)*D+CC3)*D+CC4)*D+CC5)*D+CC6)*D+CC7
      PH(D)=(((((CC8*D+CC9)*D+CC10)*D+CC11)*D+CC12)*D+CC13)*D+CC14
C      F(D)=DSQRT(POT/D)*CEXP(-CN*D+TH(-8./X))
      F(D)=DSQRT(POT/D)*CDEXP(-CN*D+TH(-8./X))
      G(D)=CDEXP(CN*D+TH(8./X))/DSQRT(TP*D)
C***
C***  FIX FOR REL. PERM.
C***
C     X=DSQRT(TPCMU*SIGL)*ROLAM
C      X=DSQRT(TPCMU*RMU*SIGL)*ROLAM
      X=SQRT(TPCMU*RMU*SIGL)*ROLAM
C***
C***  END FIX
C***
      IF (X.GT.110.) GO TO 2
      IF (X.GT.8.) GO TO 1
      Y=X/8.
      Y=Y*Y
      S=Y*Y
      BER=((((((-9.01E-6*S+1.22552E-3)*S-.08349609)*S+2.6419140)*S-32.36
     13456)*S+113.77778)*S-64.)*S+1.
      BEI=((((((1.1346E-4*S-.01103667)*S+.52185615)*S-10.567658)*S+72.81
     17777)*S-113.77778)*S+16.)*Y
      BR1=DCMPLX(BER,BEI)
      BER=(((((((-3.94E-6*S+4.5957E-4)*S-.02609253)*S+.66047849)*S-6.068
     11481)*S+14.222222)*S-4.)*Y)*X
      BEI=((((((4.609E-5*S-3.79386E-3)*S+.14677204)*S-2.3116751)*S+11.37
     17778)*S-10.666667)*S+.5)*X
      BR2=DCMPLX(BER,BEI)
      BR1=BR1/BR2
      GO TO 3
1     BR2=FJ*F(X)/PI
      BR1=G(X)+BR2
      BR2=G(X)*PH(8./X)-BR2*PH(-8./X)
      BR1=BR1/BR2
      GO TO 3
2     BR1=DCMPLX(.70710678,-.70710678)
C***
C***  FIX FOR REL. PERM.
C***
C3     ZINT=FJ*DSQRT(CMOTP/SIGL)*BR1/ROLAM
C3     ZINT=FJ*DSQRT(CMOTP*RMU/SIGL)*BR1/ROLAM
3     ZINT=FJ*SQRT(CMOTP*RMU/SIGL)*BR1/ROLAM
C***
C***  END FIX
C***
      RETURN
      END
